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ABSTRACT 

A novel method of generating artificial scintillation noise is developed and used to evaluate occupa- 
tion rates and false positive rates for surveys probing the Kuiper Belt with the method of serendipitous 
stellar occultations. A thorough examination of survey design shows that: (1) diffraction-dominated 
occultations are critically (Nyquist) sampled at a rate of 2 Fsu -1 , corresponding to 40 s -1 for objects 
at 40 AU, (2) occultation detection rates are maximized when targets are observed at solar opposition, 
(3) Main Belt Asteroids will produce occultations lightcurves identical to those of Kuiper Belt Objects 
if target stars are observed at solar elongations of: 116° < e < 125°, or 131° < e < 141°, and (4) 
genuine KBO occultations are likely to be so rare that a detection threshold of > 7 — 8cr should be 
adopted to ensure that viable candidate events can be disentangled from false positives. 
Subject headings: atmospheric effects, occultations, Kuiper Belt, Solar System: formation 



1. INTRODUCTION 

A method of detecting km-sized Kuiper Belt Ob- 
jects (KBOs) is that of serendipitous stellar occulta- 
tions (SSO). Km-sized KBOs are too faint to be ob- 
served directly, but can, in theory, be detected indi- 
rectly. The passage of a KBO through the line of sight 
to a background star will produce a characteristic per- 
turbation in the observed stellar photometric time se- 
ries, serendipitously revealing the presence of t he oc- 
culter. The method was originally proposed by Bailev 
and has been develop e d extensively jP yson 



(iy/oj, ana nas been aeveiopea extensively (Dvson 
1992J ; |Brown fe WebsteH 119971: iRoaues fe Moncuauet 



2lMlCoo rav 2003; |C^orayXSrmcr 2003; Gaud] 
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Nihei et al.l 120071). Both optical (iRoques et alj „ 
2006t IZhang et al J 12001 iBickerton et al.ll2008l hereafter 
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BKW), and x-ray (jChang et al.l 120061 : 1 Jones et al J 120061) 
occultation searches have been used to probe the popu- 
lation of the trans-Neptunian region. In this paper, we 
evaluate a variety of issues which are central to planning, 
and performing a KBO occultation survey. The methods 
we've used in our analyses are presented in Section^ and 
results for the issues examined are presented in Section[3l 
Much of the work presented herein is based on the 
planting and recovery of occultation profiles in synthetic 
photometric time series. Great care must be taken to re- 
produce the atmospheric scintillation properties of real 
photometric data, and we begin in Section 12.21 by de- 
veloping a method of generating artificial time series to 
reproduce the properties of real time series in both the 
time domain and the frequency domain. 
The SSO method was examined in the context of 
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the Kuiper Belt ([Rogues fe MoncuquetJ 12000. hereafter 
R&M) before a br eak in the KBO siz e distribution 
had been identified (Bernstein et al. 200 4]) or confir med 
(jFuentes fe HolmarJl2008t iFraser fe Kavelaarsll2009r ). In 
Section 13.11 we present updated KBO occultation rate 
estimates based on a broken power-law size distribution 
for the KBOs, and with consideration for the detectabil- 
ity of diffraction-dominated occultation events. 

Photometric time series in KBO occultation surveys 
have been obtained at a vari ety of sampling frequen- 
cies : 5 s" 1 (IZhang et all l200l). 40 a" 1 (BKW), and 45 
s- ^Roaues et al.l 12006). In Section 13.21 we show how 
Fourier methods can be used to establish a critical 
(Nyquist) sampling rate for occultation work. Numer- 
ical tests of the theory are presented in Section 13.41 

The probability an occultation will occur is propor- 
tional to the perpendicular-projected velocity of the oc- 
culter with respect to the Earth. For the Kuiper Belt, 
the motion is retrograde and the velocity is maximized 
by conducting observations at solar opposition (R&M). 
In Section 13 3| we consider the possibility that off- 
opposition observations could improve the detectability 
of smaller, and potentially more numerous objects. 

It is often useful to know the diameter of the smallest 
detectable object for data with a given signal to noise. 
In Section 13.51 we show how the signal to noise ratio, 
hereafter S/N, determines the minimum detectable size. 

Interest in small body occultations has focused largely 
on the Kuiper Belt, but other near-field objects are capa- 
ble of producing occultations. The size/distanc e degen- 
eracy for occultation events has been assessed bv lCoorayl 
(2003), with a discussion on how it may be broken for 
diffraction-dominated events. In Section 13.61 we evalu- 
ate the size/distance degeneracy further to determine the 
size/distance an occulter must have to produce an event 
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identical to a KBO occultation. We then consider the 
possibility of degeneracy with the Main Belt Asteroids 
(MB As). 

We close in Section 13.71 with an evaluation of how at- 
mospheric scintillation influences the false positive oc- 
cultation rate. We then compare false positive rates to 
the expected rates for real occultations to determine an 
appropriate detection threshold. 

Results are discussed and summarized in Sections 0] 
and [51 respectively. 

2. METHODS OF ANALYSIS 

Optimization of the observing parameters in occulta- 
tion work requires a metric by which the different pa- 
rameters can be compared. We have chosen to mea- 
sure how successfully an event can be detected with a 
'plant and recover' strategy. There are three contribut- 
ing components in this approach: a simulated occulta- 
tion lightcurve, a simulated photometric time series, and 
a detection algorithm. 

2.1. Simulating Occultation Lightcurves 

Diffraction effects become significant in an occultation 
when the occulting object has a radius smaller than the 
Fresnel scale unit, 1 Fsu = (Ad/2) 1 / 2 , where A is the 
wavelength of observation, and d is the distance to the 
occulter. At 40 AU the Fresnel scale for 550 nm light 
is 1285 m. A successful occultation lightcurve model 
for a km-sized K BO must therefore take into account 
diffraction effects dRoques et all 1 19871 iNihei et all 120071: 
iBickerton et al.l |2008[ ). We ha ve adopted a l ightcu rve 
modeling method described by iRoques et al.l (|1987f ) in 
which a circular occulting mask is assembled by placing 
progressively smaller rectangles at the periphery of the 
mask. Five "orders" of these rectangles have been used 
for all occultation lightcurve models in this work. 

All models were numerically integrated over a stellar 
disk with a projected size appropriate for the distance to, 
and spectral type of, the star being considered. Similarly, 
models were numerically integrated over a 400-700 nm 
passband. 

2.2. Simulating Photometric Time Series 

The evaluation of false positive rates ( Section |3~7|) re- 
quires that realistic artificial time series be analyzed. 
iDravins et alJ (|1997l ) showed that scintillation in a pho- 
tometric time series produces a power spectrum which 
decreases as a function of frequency. An artificial power 
spectrum can be generated and tailored to reproduce the 
properties of an observed one, and this artificial power 
spectrum can then be inverse-Fourier transformed to pro- 
duce a time series which has the same statistical proper- 
ties - both in time and frequency - as the original 'real' 
data set. 

2.2.1. The Power Spectra of Photometric Time Series 
Discrete Fourier transforms (see, for example lBracewelil 



1986) of the photometric time scries described in BKW 
were computed with a fast Fourier transform (FFT) algo- 
rithm. FFT algorithms are most efficient for time series 
composed of N = 2 n points, where n is an integer, and 
the time series were truncated to have the largest number 
of points which satisfied this requirement. The real and 



imaginary frequency components for a 409.5 s time se- 
ries (2 14 =16384 points at 40 s" 1 ) are shown in the upper 
panels of Figure [T] The amplitude of the components 
decreases as a function of frequency. Each frequency 
component was normalized by the standard deviation of 
the values within ±Q.5s , and the normalized frequency 
component distributions are presented in the lower pan- 
els of Figure [1] The frequency components are consistent 
with Gaussian variates at 95% confidence (based on re- 
duced x 2 tests). 




normalized 3te H(f) 



Fig. 1. — The noise distributions of the real and imaginary parts 
of F(w). Upper panels show the real (left) and imaginary (right) 
frequency components. The frequency components were normal- 
ized by the standard deviation of points within ±0.5s — 1 , and are 
shown as histograms in the lower left (real) and right (imaginary) 
panels with best-fit Gaussians. The frequency components are con- 
sistent with Gaussian variates based on reduced x 2 values. 

The power spectrum corresponding to Figure [1] is 
shown in Figure [2] and can be seen to obey log |F| 2 oc 
,3 log/, or \F\ 2 oc f", where \F\ is the norm of the fre- 
quency components, / is the frequency, and (3 is the 
power law slope. The power spectrum approaches a value 
of zero as the highest frequency is approached, indicat- 
ing that aliasing eff ects are likely to have been minimized 
(|Press et al.l fl992) . In observed photometric time se- 
ries, no relationship between phase and frequency was 
expected or observed. 



Fig. 2. — The power spectrum for a 40 s photometric time 
series showing a power law relationship of \F\ oc f@, with {) 1. 
The trend of the average power (0.2 s — 1 bins) is shown offset above 
for clarity. 

The photometric time series reported in BKW (56 sep- 
arate 7-25 minute segments) were Fourier transformed 
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to produce power spectra and the power law slopes 
were measured with linear least squares for frequencies: 
0.1 s _1 < / < 20 s _1 . The measured slopes were nor- 
mally distributed with (3 = —1.0, and ap = 0.3. The 
value of (3 was not observed to vary as a function of fre- 
quency. As the time series were found to contain no de- 
tectable occultation events, the structure in these power 
spectra reflects the properties of the noise. 

2.2.2. Generating Artificial Time Series 

Generation of artificial time series began with the gen- 
eration of artificial frequency components 5iiJ(/„) and 
QH(fn). The frequency components were inverse Fourier 
transformed to generate the time series. 



h k = ~— 



N-l 



N ^ 

n=0 



1 1 —2-Kikn I N 



(1) 



At a frequency, /„, the real and imaginary components 
were assigned with phases, 6 n , amplitudes, a n , and a 
normalizing constant, A, to produce a power spectrum 
with slope (5: 

5? H(f n )=f?J 2 A a n cos(0 n ), and (2) 
3 H(f n )=f^ 2 Aa n sin(e n ). (3) 

The phase angles were drawn from a uniform random 
distribution with — n < 9 < 7r, and the amplitudes, a n , 
were drawn from a Gaussian distribution with a=0, and 
<y a = 1- 

To produce a time series which is entirely real-valued, 
the property of conjugate symmetry (see, for example 
iBracewell |1986[ ) requires that the spectral components 
for the negative frequencies are the complex conjugates 
of their positive counterparts, -F(/jv-n) = F*{f n ). 

From the definition of the Fourier transform, H(0) is 
determined by the mean, h, of the time series being pro- 
duced: 



N-l 



H k =J2h n e 2 ™ kn / N , 
n=0 
JV-l 



(4) 
(5) 



n=0 



where N is the number of points. 

The normalization constant, A, was chosen to produce 
the desired noise level, er^, in the time series (see Ap- 
pendix [A| : 



A = Nah 



N-l 



-1/2 



(6) 



Frequency and time components for a 128-point arti- 
ficial data set are presented in Figure [3] to illustrate the 
method 1 . 

The variates of a true photometric time series are log- 
normally distributed; but for variances of (a^/h) 2 < 0.1, 

The code used to compute the algorithm is freely available 
from SJB. 
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Fig. 3. — Artificially-generated real (top) and imaginary (mid- 
dle) frequency components, and the corresponding time series (bot- 
tom). A complex time series is produced, but the imaginary com- 
ponents all have values of zero due to conjugate symmetry. 

the distribution is indistinguishable from a Gaussian 
(|Dravins et al.lfl997b . If the chosen noise level (ah) is 
too high, the method will produce some negative vari- 
ates, inconsistent with photometric measurements. This 
method should only be considered valid for noise levels 
of a few percent, a^jh < 0.1. 

One notable difference between observed and artificial 
photometric time series is the variability in the noise 
level. Atmospheric transparency can change on time 
scales of a few minutes during an observation, and air- 
mass changes can be present in longer observations. The 
associated changes in the flux level can be normalized by 
smoothing, but normalization of the noise level (eg. by a 
moving standard deviation) would distort the amplitude 
of an occultation signal. To model a time series contain- 
ing significant changes in the overall photometric noise 
level, the method above can be extended by modulating 
the amplitude of the artificial time series. 

If photometric noise levels are below a few percent, and 
remain stable throughout a time series, an artificial time 
series data can be generated to simulate scintillation ef- 
fects consistent with the real photometric measurements. 
The artificial time series will have the same statistical 
properties as the real data in both the time domain, and 
the frequency domain. 

2.3. The Detection Algorithm 

There are three different detection methods used 
to identify occultation event s in photometric tim e se- 
ries: the variability index (|Roques et al.l I2006D. the 
cross -correlation of template lightcurves feickcrto n et alJ 
2008), and the rank probability method (jZhang et alJ 
2008f ). The cross-correlation of template occultation 
lightcurves has been used in this work. 

Artificial time series, ij, were generated for a given 
sampling rate (Section 12. 2[) . and template occultation 
lightcurves were constructed (Section 12. ip . The occul- 
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tation lightcurves were discretely sampled at the same 
rate as the time series to form cross-correlation test ker- 
nels, ki , and the kernels were then cross-correlated with 
the time series: 



(* * k )o 



(7) 



A positive peak in (t * k)j indicates that the structure 
in the time series, t, is similar to that of the test kernel, 
k, when k is offset at position j. The detection threshold 
was set to > +8<7t*fc, the limit at which real occultation 
events are expected to occur at a higher rate than false 
positives (see Section 1377]) . 

A suite of test kernels was selected to provide coverage 
over the distances, occulter sizes, and occultation impact 
parameters of interest in the cases being examined. 

3. OCCULTATION SURVEY PARAMETERS 

Using the methods described in Section [2j we now ex- 
amine several observing parameters which are of impor- 
tance in occultation work. Numerical simulations have 
been used to determine how these parameters influence 
the rate at which detectable occultation events occur. 
Whenever possible, optimal values are established. 

We begin by determining the expected KBO occulta- 
tion rate, and proceed to examine how it is effected by 
changes in the various observing parameters. 

3.1. Expected KBO Occultation Rates 

Existing estimates of the KBO occultation rate (R&M) 
were based on a single power law size distribution, 
which has since has been succeeded by a broken 
powe r law model dKenvon fc Bromle v 2004; Pan fc Saril 
20051: iBernstein et all l2004t IFuentes fc Holmanl 120081 : 
Fraser fc Kavel aars 2009). Here, we evaluate the occul- 



tation rate expected for a KBO sky surface density based 
on the more recent broken power law model. 

An occultation at distance, d [m], with maximum de- 
tectable impact parameter, b [m], projected velocity, v 
[m s" 1 ] (projected onto the plane perpendicular to the 
line of sight), and the sky surface density, £ [deg -2 ]; will 
occur at rate: 



/' 



2^ (™ 



(8) 



The Poisson probability of observing x occultations in 
time t is: P(x) = [(fd) x /x\] e _/it , and the expected time 
between occultations is (BKW): 



ln|l -P(x > 0)| 



cxp 



(9) 



Values of 6, v, and £ are determined as follows. 

The relative perpendicula r velocity of an occ ulter in an 
uninclined circular orbit is (|Liang et al.| [20021 : 



sin 2 (e) 



+ cos(e) 
(10) 

where r^ iQ are the orbital radii of the Earth and occulter, 
e is the solar elongation of the target star, and G and 
M Q are the gravitational constant and solar mass. 



We model the sky surface density, S, with the bro- 
ken power law c umula tive size distribution described by 
iGladman et alj (|2001l . eqs. 12 & 13). Known values have 
been expressed in a single constant of proportionality, 

Qi. 

£ = Qi D k qs - 4 - 8 D 1 ^ 5 [deg" 2 ], Q 1= 3.2x10 s . (11) 

Here, Dq [km] is the diameter of the smallest KBO 
considered, and [km] is the diameter at which the 
power law breaks from the large object slope, qi, = 4.8 
i|Fraser fc Kavelaars! [200l IFuentes fc Holmanl [20081 ) . to 
the small object collisional equilibrium slope q§. 

The widths of occultation shadows (the '26' term in 
equation [9|) take different values for smaller objects in 
the diffraction dominated regime, than for la rger objects, 
for wh ich such effects can safely be ignored. iNihei et al.l 
( 2007) define shadow width by the diameter of the Airy 
ring (~ 2V3 Fsu = (6AC?) 1 ' 2 ), with a smooth transition 
between the two limiting cases: 



W = 



((6Ad)^) 



3/2 



£,3/2 



2/3 



2i?*. 



(12) 



Simulations have shown this width is a reasonable es- 
timate of the detectable width of an occultation (see for 
example, BKW, Figures 7 & 12). The cumulative surface 
density, S, includes objects larger than D, and a mean 
width, W, was used to represent the average occulter. 



2b = W 



^Xd) 1 



2 \ 3/2 ^3/2 



2/3 



where (derivation in Appendix iBj) : 



D(D>D ) = - 



(- 



dN(D) 
dD 



D dD 



2i?„ (13) 



(14) 



-is 



1-9L 



2-qs 

In 



Dn 



(ft) 



D 



(qs ^ 2) 

(qs - 2) 



The fractional exponents 2/3 and 3/2 in equation[T2"lpro- 
vide a smooth transition between terms, and the sub- 
stitution of D to obtain W gives a reliable first order 
approximation 2 . 

Times required to observe one or more occultations 
by D > Dq objects are shown in Figure [4] for opti- 
cal (A = 550nm) and X-ray (A = 0.4nm) light (the 
bands observed in occultation searches to date). The 
expected average times between events are shown (equa- 
tion [9] with P(x > 0)=0.68, or 68% confidence for one 
or more events). Curv es were computed for q a=2.0 
(jFuentes fc Holmanll2008l:lFraser fc Kavelaarj2009f ). and 
3.0. Other parameter values were: i?*=0 (point source), 
distance d=40 A U, solar elongation £= 180° (opposi- 
tion), D fc =50km (IBernstein et al] \2004 ). and q h =4.8 
(|Fraser fc Kavelaarsll2009D . 

2 As W is not a linear function of D, D cannot be directly 
substituted for D to obtain W (as in equation 113d . However, to 
a first-order approximation, W is a linear function of D, and the 
substitution is valid as an approximation. 
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Fig. 4. — Average time between occultations for KBO occultation 
searches performed in the optical (left) and X-ray (right). Different 
curves represent different possible slopes for the small object size 
distribution, qg. The waiting times required for 95% confidence 
are — In 1 0.05 ~ 3x longer (equation [9} than the average values 
shown here. 

The 95% confidence times may be more appropriate 
estimates for programs seeking to detect a single event. 
From equation^ these are — In 1 0.05 1 ~ 3x longer. 

The expected occultation rates are similar for X-ray 
and optical observations. Although the X-ray Fres- 
nel scale is small enough to offer access to the more- 
numerous lOm-sized objects, there is only a single bright 
ecliptic X-ray target (Sco X-l). The relative abundance 
of suitable optical targets (eg. M35, with 100s of upper 
main sequence stars) makes them a better choice for a 
large scale KBO occultation survey. 

3.2. Nyquist Sampling for a Diffraction-Dominated 
Occultation Event 

A measured signal should be sampled at a frequency, 
/„, which is 2x the Nyquist frequency, /n v (highest fre- 
quency to be represented). A diffraction-dominated oc- 
cultation event contains some power at all frequencies, 
and the Nyquist frequency is defined here to be the 
sampling rate for which 95% of the cumulative power 
is present in frequencies below /N y - This limit is, and 
must be, somewhat arbitrary. The goal was to identify 
an upper limit for which the overwhelming majority of 
the occultation signal strength would be represented in 
the sampled measurements. The validity of this 95%- 
power limit is supported by numerical tests which are 
presented in Section I3T41 

Power-spectra were computed for the lightcurves of 
fKBO=0-l, 0.3, and 1.0 Fsu occulters. The Fresnel-scale 
unit is a function of wavelength, A, and integration over a 
passband was simulated by averaging lightcurves within 
r KB o±15% (eg. 0.085 < r KBO < 0.115 for r KBO = 0.1 
Fsu). This tkbo ± 15% range simulates integration over 
a A±27% passband, and was chosen to represent the 400- 
700 nm passband of visible light. The lightcurves, power 
spectra, and cumulative power spectra are shown in Fig- 
ure Conveniently, the cumulative power reaches 95% 
at a wave number of k—1 Fsu -1 , regardless of tkbOi cor- 
responding to a Nyquist-sampled sampling rate of f s =2 
Fsu -1 . 

The lack of dependence on Tkbo would not generally 
be the case for tkbo "C 1 Fsu lightcurves observed in 
monochromatic light. A significant portion of the power 
would then be present in extended high-frequency ring- 
ing. When a realistic integration over a passband is con- 
sidered, the extended ringing cancels and the f^y— 1 
Fsu limit is valid, even for the tkbo =0.1 Fsu example 
shown in Figure [5] 




-2-10 1 2 
Position in Shadow (Fsu) 




0.5 1.0 

Wave Number (Fsu"' 

Fig. 5. — Diffraction-dominated occultation profiles (top) and 
the corresponding power spectra (middle), and cumulative power 
(bottom). The cumulative power reaches 95% at a frequency f^ y < 
1 Fsu -1 , and Nyquist sampling is achieved at a sampling rate 
f 3 = 2 Fsu -1 . 

For visible- light observations (A=400-700 nm) of KBOs 
at distance <i=40 AU, an occultation with a projected 
velocity of v=26 km s (ie. at solar opposition) would 
be Nyquist sampled at: 



f s = 2 Fsu" 1 • 
= 40 s _1 . 



(Ad/2) 1 / 2 



Fsu s" 1 



(15) 



Beyond 40 AU, the projected velocity asymptotically ap- 
proaches 30 km s _1 , and the sampling rate scales as 
f s oc dr 1 / 2 until the projected size of the target star 
is comparable to the Fresnel scale. 

3.3. Solar Elongation of the Target Star 

Occultations observed off-opposition have lower rela- 
tive velocities perpendicular to the line of sight, and are 
better-sampled events. But, they are less probable events 
as the occultation rate increases linearly with velocity, 
/i = 2bvT, (equation[8]) . Here, we determine how the com- 
bination of increased detectability and decreased event 
probability effect the observed occultation event rate as 
a function of solar elongation. 

Diffraction dominated occultation profiles were gener- 
ated for 50 logarithmically-spaced circular occulters with 
sizes 100 m < r K BO < 2000 m, at 40 AU, observed 
in optical light (400-700 nm). The circular occulting 
mask s were assembled out of rectangular masks (Roques 
et al. 1987; BKW), and target stars were assumed to be 
point sources. 

Artificial time series were generated with 1% 1/f noise 
(S/N = 100), as described in Section E21 Tests were 
also performed with noise levels of 2% and 4%, and gave 
similar results. Time series had 65536 points, and were 
tested with sampling rates of 10 s _1 , 20 s _1 , 40 s _1 , and 
80 s -1 (above and below the Nyquist limit described in 
Section |3~2]) . 

Relative velocities were computed at solar elongations 
between opposition and the stationary point (100 - 180° 
in increments of 5°), and the template diffraction pat- 
terns were scaled to produce templates corresponding to 
each velocity. The scaled templates were sampled at the 
rates listed above and 10 simulated events were planted 
in the time series. This was repeated for 30 evenly-spaced 
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occultation impact parameters extending to 8x the Fres- 
nel scale. The planting process is described in BKW. 

The Cross-correlation detection algorithm was used to 
attempt detection of planted events (Section 12. 3p , and 
the fraction of recovered events transitioned between ~1 
for large, easily detected events, and for events well 
below the detection threshold. The occultation rate 
was integrated by summing contributions from each Sr 
radius bin: Sfj, — 2b max (r)vS'E(r, qs) , where 6 max was 
the maximum impact parameter detected, and <5£(r, gg) 
was the sky surface density of objects in the size range 
v < tkbo < r + 5r, for a size distribution having a small- 
object slope of q$. Slopes of gs =2, and 3 were tested. The 
integral was truncated to exclude objects with tkbo > 2 
km as the sky surface density of objects larger than this 
is comparatively extremely small for the size distribution 
slopes tested. The occultation rates were normalized to 
the solar opposition value and are shown in Figure [6l Re- 
sults of the qs=3 tests were similar those for gs=2, and 
are not presented in Figure [6] 




100 110 120 130 140 150 160 170 180 
Solar Elongation (deg) 

Fig. 6. — The relative detectable occultation rate as a function 
of solar elongation (normalized to an elongation of 180°) for 10, 20, 
40, and 80 s _1 sampling rates. The reduction of D m i n (Figure 
gives some off-opposition benefit for sub-Nyquist sampling rates 
<40 s — 1 . Events are sufficiently sampled at >40 s — 1 , and any 
such advantage is lost. 

The limiting size, -D m i n , was defined to be the diameter 
at which the shadow width, 2b, reached 1 Fsu. This was 
found to be a robust measure of the diameter at which 
an object was too small to be detected. The minimum 
detectable sizes versus solar elongation are shown in Fig- 
ure [7j They were found to decrease as a function of so- 
lar elongation for the sub-Nyquist samplings (10 s _1 and 
20 s _1 ), but did not change significantly at the higher 
sampling rates. 

The improved off-opposition detectability of the more 
numerous smaller objects in the 10 s _1 and 20 s _1 tests 
suggests that detection rates in such an observing pro- 
gram would not decrease in proportion to v. Above the 
Nyquist sampling rate, the detection rate decreases in 
proportion to v, and no improvement in detectability is 
achieved by working off-opposition. 

3.4. Sampling Effects 

The Nyquist limit presented in Section 13.21 indicates 
that the theoretical optimum sampling rate is f s =2 
Fsu -1 . Here we examine the detection rate for occulta- 
tion events observed by a camera operating at different 
sampling rates. 
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Fig. 7. — The minimum detectable size as a function of solar 
elongation. The improved event sampling due to lower relative 
velocities allows detection of smaller objects for f B < 20 s — 1 , but 
has little effect at higher sampling rates. 

Occultation profiles and artificial time scries were gen- 
erated as described in Section I3~3l Our 1// noise model 
was only tested for sampling up to 40 s _1 (the sampling 
rate of our observed time series). The 1// property of 
scintillation must ul timately transition to (white) photon 
noise at > 100s _1 (jDravins et al.l[i~997r ). Our time se- 
ries were therefore simulated with Poisson variates (white 
noise) to overcome this limitation in our understanding 
of the noise beyond 40 s _1 . Though each test was per- 
formed with Poisson variates, it was repeated with the 
1// noise model for frequencies up to / = 40 s _1 , and 
results were found to be consistent with the Poisson re- 
sults for the range of frequencies tested. This simulation 
represents space-based observations very well, and is rep- 
resentative of ground-based observations up to / = 40 
s- 1 . 

Constant values of D=40 AU (ie. KBOs), and qs=2.0 
were used, and the sampling rate, /, of the camera was 
varied in logarithmic spacings 6 s _1 < / < 600 s _1 . 
The tests were bench-marked with respect to a sampling- 
specific signal to noise, S/N/. A bench mark sampling 
rate of / = 40 s _1 was used, and signal to noise ratios of 
S/N4o=100, 50, 25 were tested. At each sampling rate 
Poisson variates were drawn to represent fluxes, Ii(f), 
with mean /(/) (in photo-electrons): 

/(/)= (j)(S/N 40 ) 2 . (16) 

Target stars for occultation work are generally bright 
(V<12 mag) and sky noise was not considered. 

Simulated occultation events were planted in the data, 
and random Gaussian variates, Ri(f), with a mean, 
i?(/)=0, and standard deviation, o\r(/), were then added 
to each point in the time series to simulate read noise. 
Read noise is a function of the clamp-and-sample time 

for a pixel, t^ 2 , (|McLeadll997t ). The value, a R (f), at 
each sample rate, /, was normalized to the 40 s _1 value, 
<7 H (40): 

cT R (f) = a R (A0) (■£) ' . (17) 

It was assumed that tcs is linearly related to the full 
read-out time for the CCD. This is consistent with a 
constant-size region of a frame transfer device being read 
out at a rate which makes full use of the exposure time. 
The dead-time associated with frame transfer is typically 
1 ms and was not considered. Values of <tr(40)=0, and 
15 e~ were tested. 



KBO Occultations: Rates and Survey Design 



7 



The photon noise plus read noise signals, Ii(f)+Ri(f), 
were searched for the planted occultation events with 
the cross-correlation detection algorithm, and the small- 
est detectable size, -D m i n , was determined as described 
in Section I3~3l The density, S, of occulters larger than 
■Dmin was computed for a small object slope qs— 2 (equa- 
tion Hip , and the probability of occultation was inte- 
grated over the KBO radii as described in Section [3~3l 

The detection probabilities at different sampling rates 
are presented in Figure [5] The presence of read noise 
does not strongly influence the occultation probability 
when the photon noise is low (upper panel), but becomes 
a significant attenuating factor at higher photon noise 
levels (lower panel). This is not surprising as the read 
noise is a more significant component of the total noise 
level when photon counts are low. 
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Fig. 8. — The detectable occultation rate (normalized to peak 
value) as a function of the sampling rate. Both the influence of 
read noise (top) and that of photon noise (bottom) are shown. In 
each case, detectabil ity p eaks at the theoretical f s = 40 s" 1 rate 
predicted in Section 13.21 The presence of read noise reduces the 
detection probability at high sampling rates, and the effect becomes 
more pronounced at low signal to noise. 



3.5. Signal to Noise 

Significant increases in the occultation event rate can 
be achieved by using a larger telescope to increase the 
signal to noise in the photometric time series. Here, we 
evaluate the detectable occultation rate at different pho- 
tometric noise levels which are representative of available 
observing facilities. 

KBO (ie. 40 AU) occultation profiles were planted 
in simulated 40 s _1 (ie. optimally sampled) time scries 
containing 1/ f noise at logarithmically-spaced signal to 
noise levels of fO < S/N < fOOOO. Detection was at- 
tempted with the cross-correlation algorithm to deter- 
mine the minimum detectable size, -D min , as described in 



Section |3~!I1 Figure O shows D min versus S/N. 

The minimum diameter, in Fsu, empirically follows the 
power law: 



(S/N)- 



(18) 



with Q 2 = 3.2 ±0.2, and T] = 0.52±0.0f (by linear least- 
squares fit to logAnin = log Q2 ~ V log(5/JV) ) . TlllS 

relation was found to remain valid when tested at dif- 
ferent distances, but is expected to break down at dis- 
tances for which the target star has a projected radius >1 
Fsu. Combining the minimum detectable size with the 
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Fig. 9. — The minimum detectable size as a function of the 
photometric signal to noise. The telescope apertures required to 
achieve the S/N values presented are shown above for V=8, 10, 

ftflkAMrf e^i^cffebution for the small occulters (equa- 
tion [TT]) the sky surface density, E, can be expressed in 
terms of S/N: 



Qi Q-2 



9S D qs- 



(S/N) - r > {L - qs >. (19) 

The telescope apertures required to obtain the corre- 
sponding S/N values for apparent magnitudes of V=8, 
10, and 12 mag 3 are shown above the figure for refer- 
ence. 

3.6. Size-Distance Degeneracy and the Main Belt 
Occultation Rate 

Similar sized objects (in Fsu) at different distances will 
produce diffraction profiles that are scaled in width with 
respect to one another. This makes it possible for distant 
and nearby occulters to produce identical occultation 
lightcurves; a distant occulter will cast a wider shadow 
but the observer will move through it more quickly. 
Here, we demonstrate how such a size/distance conspir- 
acy manifests itself observationally, and the possibility 
of degeneracy with the Main Belt Asteroids (MBAs) is 
considered. 

A diffraction-dominated occultation will have a peak- 
to-peak width of (GAd) 1 / 2 « 3.5 Fsu (INihei et alJl2007L 
The duration of an occultation observed at wavelength, 
A, for an occulter at distance, d, can be expressed as: 



to 



(6Ad) 



1/2 



(20) 



where v± is the perpendicular velocity of the occulter 
described by equation[TOJ and the distance to the occulter 
is given by: 



r l + r l + 2r (B r o COS 



sin(e) 



1/2 



(21) 

Here, 7"® i0 are the Earth-Sun and occulter-Sun distances, 
and e is the solar elongation of the target star. 

Figure [TO] shows occultation time as a function of the 
distance to the occulter for targets at different solar elon- 
gations. At opposition, a diffraction-dominated occulta- 
tion by an object at 40 AU produces a shadow width 

3 Through-put of 0.65 was used to account for quantum efficiency 
(0.85), mirror reflectivity and transparency of optics (0.85), and 
atmospheric transparency (0.9). 
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similar to that observed for an occultation by an object 
at 0.08 AU, just outside of Earth's orbit. The nearby ob- 
ject would have to have the same size, in Fsu, to project 
a shadow with the same structure; and the ratio of sizes 
is given by the ratio of Fresnel scales at the two distances: 
(0.08/40) 1/2 = 0.045. Thus, a 500m diameter KBO at 
40 AU will produce an occultation shadow identical to 
that produced by a 0.045 x 500 m = 22m diameter ob- 
ject at a 0.08 AU. This is the only degenerate point for 
e = 180°; all other objects at all other distances will pro- 
duce shadows which are different in either duration or 
amplitude. 




1000 



Distance to object (AU) 

Fig. 10. — The duration of occultation as a function of the 
distance to the occulter. Different curves represent observations 
made at different solar elongations. Degeneracy of size/distance 
occurs within the Main Belt for elongations 131° < e < 141°. 

The degenerate point moves to greater distances when 
e 7^ 180. The greatest risk of degeneracy occurs when the 
degenerate point is within the asteroid belt (~ 2.0 — 3.5 
AU) . This occurs when a target is observed at elongations 
131° < £ < 141°. For elongations e < 133°, a second set 
of degenerate points enters the parameter space. These 
points represent possible prograde occultations, as the 
elongation has moved beyond the stationary point for 
the distance in question. The second set of degeneracy 
points lies within the Main Belt for elongations 116° < 
e<125°. 

At £ 136°, the degeneracy point is within the Main 
Asteroid Belt at distance d = 1.7 AU (orbital radius 2.5 
AU). There, a 500 m KBO would be degenerate with a 
(1.7AU/40AU) 1 / 2 x 500 m = 104 m MBA. 

To estimate the number of conflicting occulters we 
compute the ratio -Pmba/-Pkbo> with P = 2&wE (equa- 
tion [SJ. The shadow widths, 2b, were taken to be 3.5 
Fsu (equation [T2"]) . and velocities were computed with 
equation [10] KBO sky surface densities were calculated 
with equation [Til an d MBA sky surface densities were 
converted from the power law model of the cumulative 
luminosity function (CLF) for the Main Belt with slope 
a = 0.27 (Gladman et al. 2009, in preparation) : 

S(< to r ) = 210 x 10" ( m *- 23 ), 20 < TOR < 23. (22) 

When an object with diameter D, albedo p — 0.05, 
and heliocentric distance A is observed in reflected light 
at phase angle f3 = arcsin(dsin(e)/A) (Sine Law) at a 
distance d from the Earth, its apparent V magnitude is 
(|Coxll2000D : 



V = H{fl) + 51og(dA), (23) 

H((3)=H- 2.5 log[(l - G)*i(/3) + G$ 2 (/3)], (24) 

$ i = exp[-A i [tan( / 9/2)] B *], (25) 

H = 15.645- 2.5 \ogp- 5 log D. (26) 

The constant terms have values: Ai = 3.33, A 2 =1.87, 
£?i=0.63, and i?2 = l-22; and a slope parameter of G — 
0.2 (typical for MBAs) was used. A value of V-R=0.4 
(typical for MBAs) was then used to obtain apparent R 
magnitudes to apply in equation 1221 

The size distributions of both the KBO and MBA pop- 
ulations are modeled as power laws, and the probability 
ratio, -Pmba/-Pkbo> remains constant as a function of 
size, D, when the slopes are equal: qs = 5a + 1 = 2.35. 
In this case the probability ratio is Pmba/ -Pkbo = 2.5, 
and occultations by MBAs are expected to be more nu- 
merous. If qs < 2.35 for the Kuiper Belt, occultation 
events will be predominantly by Main Belt occulters. 

Via this transformation, a 100 m MBA has magnitude 
tor=26.6 mag, and is an extrapolation of ~3.6 magni- 
tudes beyond the observed limit of the CLF (Gladman 
et al. 2009, in preparation) The assessment remains valid 
for MBAs with sizes £>mba > 500 m, and D KB o ^ 2500 
m. 

3.6.1. Breaking the Size-Distance Degeneracy 

As the slopes for the size distributions of small KBOs 
and MBAs are not well understood, it is worthwhile to 
examine a method of breaking the size-distance degener- 
acy. 

In an occultation survey, a target star with large pro- 
jected diameter is undesirable as it causes deterioration 
of the diffraction profile. This deterioration could allow 
a shadow cast by a distant occulting object to be distin- 
guished from its near-field degenerate counterpart. De- 
terioration of a profile due to the size of the background 
star becomes significant when the star's radius is com- 
parable to the Fresnel scale. The Fresnel scale increases 
as the square-root of the distance, but the projected size 
of the background star increases linearly with distance 
to the occulter. Thus, a star which has a 'large' (~1 
Fsu) projected size at 40 AU has a relatively small (~0.3 
Fsu) projected size at 1.7 AU. The ideal target star to 
break the degeneracy will be the brightest available with 
a projected diameter of ~ 1 Fsu. 

3.7. False Positives 

The diffraction-dominated occultation shadows being 
sought in a photometric time series can be similar in 
structure to the noise, and random statistical fluctua- 
tions could be misidentified as occultations. Here we 
evaluate the false positive rate with our cross-correlation 
detection algorithm, and compare it to the expected 
occultation rate to determine an appropriate detection 
threshold. 

Artificial time series were generated with different 
noise properties (S/N, and power-spectral slope), and the 
cross-correlation detection algorithm was run on the sim- 
ulated data. Three hundred detection kernels for a point- 
source with 50m < trbo < 1200m, 10AU < c?kbo < 
160AU, and b < 3 Fsu (BKW) were used. Each time 
series had length 2 23 w 8.4 x 10 6 points (58.2 hours at 
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40 s" 1 sampling) and was guaranteed not to contain any 
real occultation events. The distribution of false posi- 
tive events in standard deviation units is presented in 
Figure [IT] with the expected Gaussian distribution for 
reference. 
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Fig. 11. — The distribution of false positives in units of standard 
deviations shown in events per second (top) , and as a fractional ex- 
cess above the rate expected for Gaussian variates (bottom). The 
use of multiple detection kernels causes the number of false posi- 
tives to be ~ 10 — 20 X higher than the rate expected for Gaussian 
variates. Changes in the signal to noise do not significantly alter 
the rates, but changing the slope of the power spectrum from -1 
(1/f noise) to (Poisson noise) increased rates by a factor of ~ 2. 



The cross-correlation detections are ~10-20x more 
common at all levels of significance because multiple de- 
tection kernels were used. But, the false positive rate 
did not increase by a factor of 300 (the number of ker- 
nels) because the kernels do not represent independent 
tests. By design, the spacing of kernel parameters tkbOj 
c^kbO) and b was chosen to allow an event to be detected 
by multiple kernels. 

The false positive rate decreases as the slope of 
the power spectrum becomes steeper. An occultation 
lightcurve has a peak in power at ~ 0.5 Fsu (see Fig- 
ured]), and a negative slope in the power spectrum indi- 
cates that more power is present in the lower frequencies. 
The relative amount of power in the noise band shared 
by occultation events decreases as the slope of the power 
spectrum steepens, and fewer false positives are observed. 

Evaluation of the false positive rate must be done on a 
per-observation basis with an accurate model of the pho- 
tometric noise, and with the same detection algorithm 
used on the real data. An excess of false positives above 
the Gaussian estimate will be produced by any method in 
which multiple passes are made over the data to extract 
different patterns. In the example presented (Figure fTT]) . 
reference to an ! n — <r' cross-correlation peak would be 
valid for the kernel tested, but would be a misleading rep- 
resentation of the statistical significance for the complete 
collection of kernels. 

3.7.1. Selection of a Detection Threshold 



There is no established detection threshold above 
which an event should be considered a viable candidate. 
We have used 8<t throughout this work and will now offer 
a justification for that limit. 

False positive events should represent a small portion 
of the number of observed events, e.g. 5% for 2a confi- 
dence. To determine a detection threshold, we consider 
the ratio of expected genuine events to false positives 
events as a function of statistical significance (in stan- 
dard deviation units). 

The false positive event rate, /if p , was computed by 
scaling the amplitude of a Gaussian to fit the false pos- 
itive distribution shown in Figure 1111 The scaling com- 
pensates for the increased false positive rate produced by 
the cross-correlation with multiple kernels. 

The expected occultation rate at each level of statisti- 
cal significance was estimated with a plant/recover pro- 
cess for a comprehensive selection of KBOs. Occultation 
profiles were tested for KBOs at 40 AU, observed in vis- 
ible light, with radii 10 m < r K BO < 2000 m (n=200, 
^?~kbo = 10 m), and impact parameters m < b < 
8000 m (n=100, 5b = 341 m). The statistical significance 
of the selected events was measured in artificial 40 s _1 
time series with S/N=25,100 and /3=1.0 (1/f noise), and 
the occultation rate for each grid point was calculated: 

5^{n, bj) = 25bv [E(n - Sr/2) - E(r< + Sr/2)] . (27) 

Each grid point was binned in statistical significance, 
and the rates, 5fi, were summed in each of the bins: 

MK) = SP ( r ^ b i)> ( a ~ <W 2 < CT k < cr + Sa/2). 

i,i 

(28) 

The occultation rates for false positives, /Zf p (<7fc), and 
expected genuine positives, /u(cfe)j are shown in FigurefT2l 
(upper panel) for different slopes of the KBO size distri- 
bution, qs=2,3. The ratio of the rates, /Lt(cfc) : /if p ((7fc), 
is shown in the lower panel. 

Variability in the expected occultation rate is due al- 
most entirely to uncertainty in the slope of the size distri- 
bution, qs (factors of > 100 x versus < 10 x for changes 
in S/N). 

From Figure [TJ] we see that the expected genuine oc- 
cultation rate exceeds the false positive rate for statistical 
significances of > la (<?s=3), or > 8c (gs=2). A limit of 
8cr has been used throughout this work. 

If the photometric noise properties are sufficiently well 
understood, the distribution of false positives can be 
modeled with artificial time series. It would then be 
possible to infer the presence of genuine events as a sta- 
tistical excess among more numerous false positives. 

4. DISCUSSION 

The parameter space for occultation studies is vast, 
including properties of the occulter (size, distance, sky 
surface density), the target star (solar elongation, pro- 
jected diameter), the time series (sampling rate, S/N, 
power spectrum slope), and the stochastic nature of the 
event itself (impact parameter). We have endeavored to 
reduce the complexity of this parameter space and de- 
termine optimum values for those parameters which are 
under an observer's control. Here, we discuss other con- 
siderations associated with SSO work. 
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Fig. 12. — The distributions of false and genuine positive occul- 
tation events in units of standard deviation (top), and the ratios 
of real-to-false positive events (bottom). False positives dominate 
for significance of < 8cr. 

The plant and recover method was adopted for a va- 
riety of tests performed, and used the same occultation 
profile for planting as was used for construction of the 
detection kernel. When used for detection in real data, 
multiple kernels are used to cover the parameter space of 
interest. The parameters of an occultation event (tkbo, 
^kbOj and b) would not be expected to match any of the 
kernels exactly, but would be recovered due to overlap 
in their ranges of sensitivity. Use of the same profile for 
both planting and detection improved the efficiency of 
our testing code, but does not compromise the validity 
of any results. 

Our finding that occultation events are best sampled at 
/ s =2 Fsu" 1 is based on 95% of the cumulative power be- 
ing in frequencies below /n y =1 Fsu - . This is true only 
for diffraction profiles integrated over a range of wave- 
lengths. An occultation profile observed in monochro- 
matic light will contain extended ringing and retain con- 
siderable power at frequencies beyond l Fsu -1 . The sam- 
pling rate we provide here is valid for filters with ~ 100 
nm pass-bands, typical of broad-band filters used in op- 
tical astronomy. 

Estimates of the occultation rate as a function of the 
sampling rate were based on plant /recover tests in time 
series of Poisson variates rather than 1// noise. Our pho- 
tometric time series were sampled at 40 s _1 and we did 
not extrapolate the power spectrum to higher frequen- 
cies in order to simulate the noise at higher sampling 
rates. The tests were repeated with 1// time series for 
frequencies up to /=40 s _1 , and these were found to be 
consistent with the Poisson noise tests for the range of 
frequencies compared. 

The detectable occultation rates shown in Figure [4] are 
considerably lower than earlier estimates. An occultation 
by a D «1000 m object would be expected at a rate 
of /i = 10~ 8 — 10~ 10 s _1 , where the previous rate was 
5 x 10 -7 s _1 (R&M). The rates presented here take into 
account a break in the size distribution which was not 
well established when earlier estimates were made. The 



break to a shallower size distribution at tkbo ~ 25 km 
dramatically reduces the expected sky surface density of 
km-sized KBOs. 

The cross-correlation detection algorithm was used in 
an occultation search of 40 s _1 photometric time series 
from two B9V stars in M35 (BKW) . In real photometric 
time series, false positives were detected at a higher rate 
than would be predicted in Section [3~7l of this work. The 
critical difference between observed and simulated time 
series is that photometric stability cannot be maintained 
for an observed time series. Atmospheric transparency 
can change on time scales of a few minutes, and thereby 
change the signal to noise ratio of the photometry. To 
guard against this, detection algorithms should normal- 
ize a candidate event's statistical significance by a more 
local measure of the standard deviation than has been 
used in the past - within seconds, rather than minutes 
of a candidate point. 

5. SUMMARY 

We have provided a comprehensive assessment of the 
various observing parameters relevant in an SSO survey. 
Our aim has been to identify optimal values for parame- 
ters to guide future observations in this field. The results 
we have obtained are summarized here. 

• We have estimated occultation event rates based 
on a broken power law size distribution, and found 
rates of 10~ 8 — 10~ 10 s _1 for km-sized objects. The 
variability is due to the large uncertainty in the size 
distribution slope, qs, for small objects. 

• Fourier analysis indicates that 95% of the cumula- 
tive power in the power spectrum of a A = 400—700 
nm occultation profile is present in frequencies of 
/ < 1 Fsu _1 , indicating that an event can be crit- 
ically (Nyquist) sampled at a rate of f s =2 Fsu -1 . 
To probe the Kuiper Belt at 40 AU, this corre- 
sponds to a sampling rate of 40 s" 1 . Numerical 
tests of the sampling rate indicate that sampling 
significantly above the Nyquist limit (>80-100 s _1 ) 
increases read-noise, and reduces the detectable 
KBO occultation rate. 

• The photometric noise present in 40 s _1 ground- 
based observations has a power spectrum with the 
form 1/f 13 , with /3 » 1. We have used Fourier 
methods to generate 1// noise artificially and sim- 
ulate realistic photometric time series. Plant-and- 
recover tests were performed with these time series 
to numerically test various modes of observation. 

• When optimally sampling, the detectable occulta- 
tion event rate is maximized by observing a target 
at solar opposition (within ±20 — 30°). Observ- 
ing off-opposition does not significantly improve de- 
tectability of smaller, more numerous objects un- 
less the sampling rate of the camera is sub-critical 
for an observation at opposition. 

• In a diffraction-dominated occultation event which 
is well sampled, the occulter distance will manifests 
itself in the width (duration) of the event. A given 
distance is degenerate with either one, or three oth- 
ers, depending on the solar elongation of the obser- 
vation. Occultation events by KBOs at 40 AU are 



KBO Occupations: Rates and Survey Design 



11 



degenerate with 2.0-3.5 AU MB As for observations 
performed at elongations of 116° < e < 125°, or 
131°<e<141°. 

• Occultations by MBAs may occur at a rate compa- 
rable to those for KBOs depending on the slopes of 
the size distributions for the two populations. An 
MBA occultation observed at solar opposition will 
be ~ 0.5 x the duration of a KBO event, and the 
two should be distinguishable. The sampling rate 
adopted for a survey targeting multiple distances 
should be the highest of the critical sampling rates 
for each individual distance (~ 80 s~ x for MBAs). 

• Due to the expected rarity of KBO occultation 
events, false positives are likely to be more common 
than real events for levels of statistical significance 
< 7 — 8a, depending on the slope of the KBO size 
distribution. 
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APPENDIX 

DERIVATION OF THE 1/F-NOISE NORMALIZATION CONSTANT 

The normalization constant, A (equations |2J and [3]) is need to produce an artificial time series with 1/f scintillation- 
like properties. With variables as defined by equations [2] through \5\ in Section [2.2.21 the value of A can be determined 
as follows: _ 

Taking the time variate hk = h + £k (e = 0, a e = ah), substitution of equation [5] into Parseval's theorem yields: 
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THE AVERAGE OBJECT DIAMETER IN A BROKEN POWER LAW SIZE DISTRIBUTION 

Equation [T3] describes the average diameter of an object in a broken power law size distribution. The expression can 
be derived as follows: 
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where N is as defined in equation ??. These two integrals can be expressed as: 
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With substituting of equations IB2i and IB3I equation IB II can be reduced to give: 
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Equation IB4I is undefined for gL = 2, but this is well below current estimates of its value (gL ~ 4.8 lFuentes Sz Holmanl 
2008; Frascr & Kavelaars 2009). The equation, as shown, is also undefined for qs = 2, but has a finite limiting value. 
It can then be rearranged to isolate the (gs — 2) terms, and take the limit as gs — > 2: 
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The limit can then be computed with L'Hopital's rule by taking derivatives of the numerator, (1 — (Do/D k ) qs 2 ), 
and denominator, (2 — qs), in the second term of Equation IB5I 
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